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Abstract 



Finite element analysis has been used successfully to estimate the effective properties of many types of composites. 
The prediction of effective elastic moduli of polymer-bonded explosives provides a new challenge. These particulate 
composites contain extremely high volume fractions of explosive particles (> 0.90). At room temperature and higher, 
the Young’s modulus of the particles can be 20,000 times that of the binder. Under these conditions, rigorous bounds 
and analytical approximations for effective elastic properties predict values that are orders of magnitude different 
from the experimental values. In this work, an approach is presented that can be used to predict three-dimensional 
effective elastic moduli from two-dimensional finite element simulations. The approach is validated by comparison 
with differential effective medium estimates and three-dimensional finite element simulations. The two-dimensional 
finite element approach has been used to determine the properties of models of polymer bonded explosives and PBX 
9501 in particular, containing high volume fractions of circular and square particles with high modulus contrasts. 
Results show that estimates of effective elastic properties from two-dimensional finite element calculations are close 
to the values predicted by the differential effective medium approach for a large range of volume fractions and modulus 
contrasts. Two- and three-dimensional finite element estimates for volume fractions from 0.70 to 0.90 and found not to 
differ considerably. Simulations of models of polymer bonded explosives and PBX 9501 show that the microstructure, 
the amount of discretization, and the type of element used play a considerable role in determining the value predicted 
by finite element simulations. The effective elastic moduli of PBX 9501 predicted by finite element calculations can 
vary from 200 MPa to 10,000 MPa depending on the microstructure and level of discretization used. The results also 
suggest that if a microstructure can be found that accurately predicts the elastic properties of PBX 9501 at a certain 
temperature and strain rate, then the same microstructure can be used to predict elastic properties at other temperatures 
and strain rates. 

Keywords : particle-reinforced composites, microstructure, modeling, elastic properties 

1 Introduction 

Experimental determination of the mechanical properties of polymer-bonded explosives (PBXs) is hazardous and 
highly expensive. An alternative to experimentation is the determination of material properties using micromechanics 



2 




methods. However, rigorous bounds and analytical approaches predict values of elastic modulus that are orders of 
magnitude different from the experimental values. Hence, detailed numerical simulations are required for the predic- 
tion of effective properties of PBXs. The goal of this work is to investigate if the elastic properties of PBXs can be 
determined using finite element analyses of PBX microstructures. 

PBXs are particulate composites containing explosive particles and a rubbery binder. The particle size distributions 
in these materials and the mechanical properties of the constituents are required prior to finite element simulations 
of the microstructure of these materials. The size distribution of the explosive particles in PBXs can be obtained 
experimentally. The mechanical properties of the particles can be determined from molecular dynamics simulations. 
Mechanical testing of the rubbery binder does not involve any risk of explosion. Hence, mechanical properties of the 
binder can be determined from laboratory tests. 

The major challenges involved in finite element modeling of PBXs are the high volume fraction of particles in 
PBXs ( f p > 0.90) and the high modulus contrast between particles and binder (E p /Et, = 10,000 to 20,000) at and 
above room temperature and at low strain rates. Because of the high volume fractions, microstructures are difficult 
to generate and even more difficult to digitize from micrographs. In addition, the high modulus contrast leads to 
large stress concentrations at the interface of particles and the binder. These stress concentrations cannot be simulated 
adequately without a high degree of mesh discretization. Hence, three-dimensional finite element models of these 
materials require considerable computational power. 

In this work, a procedure is outlined whereby two-dimensional finite element simulations can be used to estimate 
the three-dimensional elastic moduli. In order to verify the accuracy of the approach, finite element estimates on unit 
cells containing circular particles are compared with differential effective medium estimates [1, 2] of effective elastic 
moduli for a wide range of volume fractions and modulus contrasts. Two-dimensional finite element estimates for 
selected microstructures are then compared to three-dimensional estimates for particle volume fractions of 0.70, 0.75 
and 0.80. 

Next, an appropriate model of PBX 9501 is determined from simulations of manually generated and randomly 
generated microstructures based on the particle size distribution of PBX 9501. Finite element calculations are then 
performed on this model to estimate the elastic moduli of PBX 9501 for various temperatures and strain rates. Finally, 
the estimated moduli are compared with experimental data to gauge the applicability of the approach to polymer 
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bonded explosives. 



2 Finite element approach 

Numerical determination of the effective properties of particulate composites involves the calculation of the stress and 
strain fields for a representative volume element (RVE) that simulates the microstructure of the composite. These 
stresses and strains are averaged over the volume (V) of the RVE. The effective elastic stiffness tensor Cg k[ of the 
composite is calculated from the tensor relation 

[ UijdV = qf w [ e kl dV, (1) 

Jv Jv 

where are the stresses and tij are the strains. 

The stress-strain relation for the RVE can be written in Voigt notation as 
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where { a) v is the volume average of the quantity a. Equation (2) can also be inverted and written in terms of the 
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directional Young’s moduli, Poisson’s ratios and shear moduli as 
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where the directional Young’s moduli Eft, f, and E % f differ by a small amount from the isotropic Young’s modulus 
-E e ff- Similarly, the mean of the Poisson’s ratios v[ f, vf[, i/ff , zzgf , , and is assumed to be close to the isotropic 

Poisson’s ratio p e ff and the mean of the shear moduli G\ f , G|f , and G'fj is assumed to be close to the isotropic shear 
modulus G e ff. The bulk mechanical response of a random particulate composite is isotropic. However, the size of a 
RVE of such a composite that can be simulated numerically is necessarily much smaller than the bulk and may not be 
isotropic. In this work, it is assumed that the deviation from isotropy of an RVE is small. 

The determination of effective isotropic elastic modulus of particulate composites using three-dimensional fi- 
nite element analysis therefore requires the determination of volume averaged stresses and strains under appropriate 
boundary conditions. Once these stresses and strains are known for an RVE, equation (3) can be used to determine 
the effective elastic moduli of the RVE. An investigation into effects of boundary conditions revealed a uniaxial state 
of stress develops in a RVE when a uniform boundary displacement is applied in one direction and periodicity of 
boundary displacements is maintained. This means that the effective Young’s modulus and Poisson’s ratio in a partic- 
ular direction can be directly computed from equation (3) by applying a uniform displacement in that direction. The 
effective isotropic Young’s modulus and Poisson’s ratio can then be calculated by averaging the values obtained in the 
three orthogonal directions from three simulations. The effective shear modulus can be calculated directly from the 
effective isotropic Young’s modulus and Poisson’s ratio. 

The estimation of effective elastic properties from three-dimensional simulations is straightforward in principle. 
However, not only is it difficult to generate RVEs containing high volume fractions of random spherical particles, 
meshing such three-dimensional RVEs is nontrivial. Even if RVEs can be generated and meshed, solution of the 
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finite element system of equations can be prohibitive in terms of computational cost. A two-dimensional approach is 
described below that requires relatively little computational expense. This approach can be used to arrive at reasonable 
estimates of three-dimensional effective properties of particulate composites. 

In the two-dimensional simulations of particulate composites performed in this study, a cross-section of the com- 
posite is chosen that contains the same area fraction of particles as the volume fraction of particles in the composite. 
Finite element analyses are then performed on this square RVE assuming a state of plane strain. The stress-strain 
equation (3) can be written in planar or two-dimensional form as 
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where Ef[, E% f, i> 2 i, and i'j'IJ are two-dimensional Young’s moduli and Poisson’s ratios and different from their 
three-dimensional counterparts. A brief explanation of two-dimensional elastic moduli is given in Appendix A. 

Uniform displacements are applied to the ‘ 1 ’ and ‘2’ directions of the RVE and periodic displacements are enforced 
along the remainder of the boundary to arrive at uniaxial average stresses and the corresponding average strains in the 
RVE. The effective moduli are calculated using equation (4) and the directional moduli are averaged to arrive at a 
mean effective Young’s modulus (E^f) and a mean effective Poisson’s ratio (v^g) for the RVE. These moduli can 
also be thought of as effective transverse moduli of a composite containing unidirectional cylinders. 

It is assumed that the transverse moduli are equal to the two-dimensional apparent moduli of an isotropic, homoge- 
neous material subjected to plane strain. Since equation (4) is a two-dimensional stress-strain law, the two-dimensional 
moduli have to be converted to three-dimensional moduli using the relations [3] 



^eff = V&IQ- + V&) 

E e g = i?eff [1 - (t^eff) 2 ], (5) 



where u e g is the effective Poisson’s ratio in three dimensions and E e g is the three-dimensional effective Young’s 
modulus. 
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3 Validation of approach 



3.1 Two-dimensional FEM vs. differential effective medium estimates 

The differential effective medium (DEM) approach [1, 2] has been found to provide excellent estimates of elastic 
properties of particulate composites [2]. To determine if the two-dimensional finite element approach provides rea- 
sonable estimates of three-dimensional effective elastic properties, finite element analyses were performed on the 
two-dimensional RVEs containing 10% to 92% by volume of circular particles that are shown in Figure 1. The par- 
ticles were assigned a Young’s modulus of 100,000 MPa and a Poisson’s ratio of 0.2. The Young’s modulus of the 
binder was varied from 1 MPa to 10,000 MPa in multiples of 10 and the binder Poisson’s ratio was chosen to be 0.49. 
The effective properties predicted by the finite element approach were then compared with DEM estimates. 

In the RVEs shown in Figure 1, circular particles were placed sequentially at random locations in each RVE such 
that no two particles were in contact. The particle size distribution in each RVE was based on that of the dry blend of 
PBX 9501 [4], The planar moduli obtained from the finite element calculations were converted into three-dimensional 
moduli using equations 5. Figure 2(a) shows effective three-dimensional Young’s moduli of the RVEs from finite 
element (FEM) and differential effective medium (DEM) calculations. Figure 2(b) shows the corresponding three- 
dimensional Poisson’s ratios. The solid lines in the figures are lines of constant modulus contrast between particles 
and binder and represent the effective elastic property calculated using DEM. The finite element results are shown as 
squares. 

The FEM and the DEM predictions are seen to agree remarkably well for all modulus contrasts and for all the 
simulated volume fractions of particles, implying that the two-dimensional FEM approach described in the previous 
section gives reasonable estimates of effective elastic properties over a considerable range of volume fractions and 
modulus contrasts. 

3.2 Two-dimensional vs. three-dimensional FEM estimates 

To obtain a more direct estimate of the effectiveness of the two-dimensional FEM approach, a second set of finite 
element simulations was performed on two- and three-dimensional RVEs containing particle volume fractions of 0.70, 
0.75, and 0.80. The two-dimensional RVEs are shown in Figure 3(a) and the three-dimensional RVEs are shown in 
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Figure 3(b). The Young’s modulus used for the particles was 15,300 MPa and the Poisson’s ratio was 0.32. The 
Young’s modulus of the binder was 0.7 MPa and the Poisson’s ratio was 0.49. These values correspond to the moduli 
of HMX and the binder of PBX 9501 at room temperature and low strain rate [5, 6], 

For the two-dimensional simulations, each RVE was discretized using six-noded triangles and approximately 
60,000 nodes. Displacement boundary conditions were applied and the average stresses and strain were used to 
calculate two-dimensional effective Young’s moduli and Poisson’s ratios. The two-dimensional moduli were then 
converted into three-dimensional moduli using equations 5. The three-dimensional RVEs were discretized using 10- 
noded tetrahedral elements with a total of around 100,000 nodes. A uniform displacement was applied perpendicular 
to one of the faces of the RVE while the remaining faces were constrained to displace in a periodic manner. The 
three-dimensional Young’s modulus and Poisson’s ratio were directly calculated from the average stresses and strains 
for these three-dimensional models. 

Figure 4(a) shows the predicted effective Young’s moduli of the two- and three-dimensional RVEs shown in Fig- 
ure 3. Though the three-dimensional models predict higher values of Young’s modulus than the two-dimensional 
models, the difference between the two- and three-dimensional estimates of Young’s modulus is small relative to the 
modulus contrast between the particles and the binder. Figure 4(b) shows the effective Poisson’s ratio for the two- 
and three-dimensional finite element models. The three-dimensional models predict lower values of Poisson’s ratio 
than the two-dimensional models. Both two- and three-dimensional models predict a sharp decrease in Poisson’s ratio 
with increasing particle volume fraction. These results show that the two-dimensional FEM approach can be used to 
obtain acceptable estimates of effective elastic properties of random particulate composites and hence as a substitute 
for detailed three-dimensional calculations. 



4 Modeling polymer bonded explosives 

A micrograph of the polymer bonded explosive PBX 9501 [7] is shown in Figure 5. The HMX particles are shown to 
be irregularly shaped and present in a large number of sizes. Two length scales can be identified from the micrograph. 
The first is the scale of the larger particles that occupy most of the volume. The second scale is that of the particles 
filling the interstitial spaces between the larger particles. Because of these two different length scales, it is extremely 




difficult to use a digital image at a single scale to generate two-dimensional microstructures of PBX 9501 or other 
PBXs. 

Most micromechanical calculations for PBXs have been carried out using subgrid models that use simplified 
models of the microstructure (for example, spherical grains coated with binder or spherical voids in an effective PBX 
material [8]). Closed form solutions from these simple models have been used to provide properties for macroscopic 
simulations. More detailed calculations have used microstructures containing ordered arrays of circles or polygons in 
two or three dimensions to model PBXs [9, 10]. These models do not reflect the microstructure of PBXs and hence 
have limited use for predicting thermoelastic properties. Better two-dimensional approximations of the microstructure 
have been constructed from digital images of the material and used by Benson and Conley [11] to study some aspects 
of the micromechanics of PBXs. However, such microstructures are extremely difficult to generate and require state- 
of-the-art image processing techniques and excellent images to accurately capture details of the geometry of PBXs. 
More recently Baer [12] has used a combination of Monte Carlo and molecular dynamics techniques to generate three- 
dimensional microstructures containing spheres and oriented cubes that appear to represent PBX microstructures well. 
However, the generation of a single realization of these microstructures is very time consuming and often leads to a 
maximum packing fraction of about 70%. Periodicity is also extremely difficult to maintain in the RVEs generated by 
this method. 

One of the goals of this work was to explore a fast, yet accurate, means of predicting the effective elastic properties 
of PBXs. With this goal in mind, accurate representation of the microstructure of PBXs was foregone. Instead, the 
explosive crystals were modeled using circles or squares. 

Both manually and automatically generated two-dimensional RVEs of PBX 9501 are discussed in this section. 
The effective elastic moduli of these RVEs from finite element calculations are then compared with experimentally 
determined moduli of PBX 9501. As in the previous section, the particles in the RVEs represent HMX and have a 
Young’s modulus of 15,300 MPa and a Poisson’s ratio of 0.32 [5] and the binder has a Young’s modulus of 0.7 MPa 
and a Poisson’s ratio of 0.49 [6]. The effective properties computed from these component properties are compared 
with the Young’s modulus of 1,013 MPa and Poisson’s ratio of 0.35 [6] of PBX 9501. 
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4.1 Manually generated microstructures 



Figure 6 shows six two-dimensional microstructures representing PBX 9501. Each RVE is filled with particles of 
various sizes since PBXs are typically a mixture of coarse and fine grains with the finer grains forming a filler between 
coarser grains. The volume fractions of particles in each of these models is 90±0.5%. All the models possess square 
symmetry. The RVEs have also been designed so that particle-particle contact is avoided. Finite element simulations 
were performed on these RVEs using six-noded triangular elements such that the circular particles were discretized 
accurately. The approach discussed in Section 2 was used to determine the effective Young’s modulus and Poisson’s 
ratios of these RVEs. 

Table 1 shows the effective Young’s modulus and Poisson’s ratios of the six RVEs. The values of the effective 
Young’s modulus range from 42 MPa to 192 MPa, with a mean of 132 MPa and a standard deviation of 54 MPa. In 
comparison, the Young’s modulus of PBX 9501 is 1,013 MPa; around 500% to 800% of that predicted using the six 
RVEs. Among the six RVEs, the Young’s modulus appears to be higher for the models with lower amounts of binder 
along the edges of the RVE. Models 1 through 3 have a single large particle and many smaller particles and show 
approximately the same effective behavior. Model 4, with a smaller ratio between the radius of the largest and the 
smallest particles, has a lower effective Young’s modulus than the mean value. 

The effective Poisson’s ratios of the six RVEs range from 0.25 to 0.34, with a mean of 0.33. The mean Poisson’s 
ratio of the RVEs is not much different from the measured value for PBX 9501. For models 1, 2, and 3, the effective 
Poisson’s ratio is within 10% of that of PBX 9501. However, the Poisson’s ratio of model 4 is 25% higher and those 
of models 5 and 6 are 20% and 28% lower than the Poisson’s ratio of PBX 9501, respectively. In comparison, the 
two-dimensional Poisson’s ratios of models 4, 5, and 6 are 0.79, 0.39, and 0.33, respectively, whereas that of PBX 
9501 is 0.54. Since these two-dimensional values reflect the actual strains obtained from the two-dimensional finite 
element calculations, it is clear that when compared to PBX 9501, model 4 is too compliant while models 5 and 6 are 
too stiff. These observations are also attributable to the amount of binder along the edge of a RVE. 

On average, the FEM estimate of Young’s modulus is around 13% of the experimentally determined value for PBX 
9501. The reason for this difference could be that preferential stress paths (stress-bridging) in the microstructure of 
PBX 9501 lead to increased stiffness. Results shown in Figure 4 indicate that the difference between estimates from 
two- and three-dimensional models is not large. Hence, it is unlikely that the considerably lower stiffness of the six 
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model RVEs compared to PBX 9501 is due to the use of two-dimensional models. 

It is possible that the Young’s modulus is underestimated by the six RVEs because the volume fraction of particles 
in each of these microstructures is lower than 92%. Figure 7(a) shows a model containing 92% particles by volume 
that is a modified version of model 6 in Figure 6. The effective Young’s modulus of the RVE containing 92% particles 
is 218 MPa and the effective Poisson’s ratio is 0.28. In comparison, model 6 contains around 90% particles by volume 
and has an estimated Young’s modulus of 192 MPa and a Poisson’s ratio of 0.25. Though there is some increase in the 
Young’s modulus due to increase in particle volume fraction, the value is still around 20% that of PBX 9501. Hence, 
the volume fraction cannot be the only factor leading to underestimation of the Young’s modulus by the manually 
generated microstructures. 

If a 256x256 square grid is overlaid on top of the model shown in Figure 7(a) and cells in the grid are assigned 
HMX properties if they contain more than 50% particle and binder properties otherwise (referred to as the ’’square 
grid overlay method”), the new microstructure shown in Figure 7(d) is obtained. The RVE in Figure 7(d) was modeled 
using 256x256 four-noded square finite elements. The modified model predicts a Young’s modulus of 800 MPa and a 
Poisson’s ratio of 0.14. The modulus is now closer to that of PBX 9501 but the Poisson’s ratio is much lower. 

The stiffer response of the model shown in Figure 7(d) is partially due to increased stiffness of four-noded finite 
elements. To eliminate the possibility that the stiffer response of the modified model is due to element locking caused 
by the nearly incompressible binder, the binder was also modeled as a two-parameter Mooney-Rivlin rubber instead of 
a linear elastic material. An incremental analysis using the Mooney-Rivlin material model yielded values of Young’s 
modulus and Poisson’s ratio that were within 5% of those for the linear elastic models. This result implies that the 
difference between the effective properties of the RVEs in Figures 7(c) and 7(d) is primarily due to the slight difference 
in microstructure. The same behavior has also been observed for all the models shown in Figure 6. 

The model in Figure 7(d) has a Young’s modulus that is four times higher than that of the model in Figure 7(c). This 
result suggests that the square grid overlay method can lead to the automatic incorporation of preferential stress paths 
into RVEs containing high volume fractions of particles. In addition, it is easier to generate a square mesh and use the 
square grid overlay method than to model circular particles using triangular elements. Hence, the square grid overlay 
method is used to discretize the randomly generated microstructures based on the actual particle size distribution of 
PBX 9501 discussed in the next section. 
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The manually generated microstructures were based on a qualitative assessment of the particle size distribution of 
PBX 9501 without regard for the actual size distribution. In addition, particles in the RVEs were not allowed to come 
into direct contact leading to a gross underestimation of the effective Young’s modulus. The first of these issues is 
addressed in the following section where the actual particle size distribution of PBX 9501 is used to generate model 
microstructures. The particle contact problem is addressed by approximating the model microstructures using the 
square grid overlay method. 

4.2 Randomly generated microstructures 

The preferred method for generating close packed microstructures from a set of particles is to use Monte Carlo based 
molecular dynamics techniques [13] or Newtonian motion based techniques [14]. In both methods, a distribution of 
particles is allocated to the grid points of a rectangular lattice using a random placement method. Molecular dynamics 
simulations or Newtonian dynamics calculations are then carried out on the system of particles to reach the packing 
fraction that corresponds to equilibrium. A weighted Voronoi tessellation [15] is then carried out on the particles with 
the weights determined by the sizes of the particles. The particles are next moved towards the center of the packing 
volume while maintaining that they remain inside their respective Voronoi cells. The process is repeated until all the 
particles are as tightly packed as possible. Periodicity of the particles at the boundaries is maintained by specifying 
extra particles at the boundaries that move in and out of the volume. This process, with some modifications, has been 
the only efficient method of generating close packed systems of particles in three dimensions. However, it is difficult to 
get packing fractions of more than 70-75% when using spherical particles. It is also quite time consuming to generate 
a tight packing. 

In two dimensions, a faster method can be used for generating particle distributions, that is random sequential 
packing. The largest particles are placed randomly in the volume followed by progressively smaller particles. If there 
is any overlap between a new particle and the existing set, the new particle is moved to a new position. If a particle 
cannot be placed in the volume after a certain number of iterations, the next lower sized particle is chosen and the 
process is continued until the required volume fraction is achieved. Though this method does not preserve the particle 
size distributions as accurately as the Monte Carlo based molecular dynamics methods, it is much faster and can be 
used to generate high packing fractions in two dimensions without particle locking. In three dimensions, this method 
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is highly inefficient and packing above 60% is extremely time consuming to achieve. 



4.2.1 Circular particles - PBX 9501 dry blend 

HMX particle size distributions of the dry blend of PBX 9501 have been listed by Wetzel [6] and Rae et al. [16]. The 
coarse and fine particles are blended in a ratio of 1 :3 by weight and compacted to form PBX 9501 . Four microstructures 
based on the particle size distribution of the dry blend are shown in Figure 8. The number of particles used for the 
four microstructures are 100, 200, 300, and 400. The RVE widths are 0.65 mm, 0.94 mm, 1.13 mm and 1.325 mm, 
respectively. The particles occupy a volume fraction of about 85-86%. The remaining volume is assumed to be 
occupied by a mixture of binder and fine particles of HMX that are well separated in size from the smallest particles 
shown in the RVEs. This ‘dirty’ binder contains around 36% particles by volume. 

Particles in the RVEs were assigned HMX properties and the elastic properties of the dirty binder were calculated 
using the differential effective medium approximation [2], The Young’s modulus of the dirty binder used in the 
calculations was 2.0816 MPa and the Poisson’s ratio was 0.4813. Two sizes of square grids were used for overlaying 
on the RVEs: 256x256 and 350x350. The square grid overlay method was applied to assign materials to grid cells. 
Each cell in the grid was modeled using a four-noded finite element. Table 2 shows the effective Young’s modulus and 
Poisson’s ratio predicted for the four models of the dry blend of PBX 9501. 

The FEM calculations overestimate the effective Young’s moduli of PBX 9501 by 200% to 400% if the 256x256 
grid is used. If the 350x350 grid is used the effective Young’s moduli vary from 100% to 300% of that of PBX 
9501. The effective Poisson’s ratio is underestimated considerably in both cases. The stiffer response of these RVEs 
is partly due to the creation of continuous stress paths across a RVE when a regular grid is overlaid on the RVE and 
partly because four-noded elements are inherently stiffer [17]. It can also be seen that the effective Young’s modulus 
increases as the number of particles in the RVE increases. This is because the same grid has been used for the four 
microstructures, leading to poorer resolution of the geometry with increasing particle density per grid cell. 

4.2.2 Circular particles - pressed PBX 9501 

Figure 9 shows four RVEs based on the particle size distribution of pressed PBX 9501 [7], The pressing process leads 
to particle breakage and hence a larger volume fraction of smaller sized particles. Fewer larger sized particles remain 
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as is reflected in the generated microstructures containing 100, 200, 500, and 1000 particles. In this case, the RVE 
widths are 0.36 mm, 0.42 mm, 0.535 mm, and 0.68 mm, respectively. Thus, the 1000 particle RVE for the pressed 
piece has dimensions similar to the 100 particle RVE for the dry blend. The size of the RVE that can be adequately 
discretized is therefore smaller for the pressed piece. Each RVE was discretized into both 256x256 and 350x350 
four-noded elements. Each element was assigned material properties according to the square grid overlay method. 

The particles in the RVEs occupy volume fractions from 0.86 to 0.89 and the target volume fraction of 0.92 is not 
attained in any of the RVEs. A dirty binder, whose properties were calculated using the differential effective medium 
approach, was used in the effective stiffness calculations. The effective moduli of the dirty binder for the four models 
are shown in Table 3. The effective Young’s modulus and Poisson’s ratio of the four RVEs from FEM calculations are 
shown in Table 4. The 256 x 256 element models overestimate the Young’s modulus of PBX 9501 by factors increasing 
from 3 to 6 with increasing RVE size. The estimates from the 350 x 350 element models are lower but still 2 to 5 times 
higher than the Young’s modulus of PBX 9501. The estimated Poisson’s ratios are quite low compared to that of PBX 
9501. 

The four models based on pressed PBX 9501 predict Young’s moduli that are 1.5 to 2 times higher than values 
predicted by models based on the dry blend of PBX 9501 . For the 100 and 200 particle pressed PBX 9501 models, the 
single large particle contributes considerably to the stiffer response. For the 500 and 1000 particle pressed PBX 9501 
models, errors in the discretization of particle boundaries lead to additional stress bridging paths and hence a stiffer 
response is obtained. 

4.2.3 Square particles - pressed PBX 9501 

With an increase in particle volume fraction, there is an increase in the number of particle -particle contacts in par- 
ticulate composites. If the particles in a RVE of such a material are circular, triangular elements cannot be used to 
discretize the RVE because of poor element shapes in regions of contact. Additionally, a rectangular grid cannot rep- 
resent the geometry of circular particles accurately. An alternative is to use square particles that are aligned with a 
rectangular grid to represent the microstructure. 

The particles shown in the three microstructures in Figure 10 are based on the size distribution of pressed PBX 
9501. These distributions have been generated by placing square particles in a random sequential manner in a square 
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grid. The smallest particles occupy a single subcell of the grid. Larger particles are chosen from the particle size 
distribution so that they fit into an integer multiple of the grid size. 

In the three models shown in Figure 10, the particle size distribution for the pressed piece is truncated so that the 
smallest (< 30 microns) particle sizes in the distribution are not used in order that the grid size does not become too 
large. The RVEs are filled with particles to volume fractions of about 86-87%. The remaining volume is assumed 
to be occupied by a dirty binder. The number of particles in the first model is 700 and the RVE width is around 3.6 
mm. The second model contains 2,800 particles and the RVE width is about 5.3 mm, while the third model contains 
1 1 ,600 particles and has a width of 9 mm. Note that these RVEs are considerably larger than those used for the circular 
particles. Smaller RVEs are not used because of the difficulties associated with fitting particles into integer multiples 
of subcell widths. The properties of the dirty binder for the three models are shown in Table 5. 

Finite element calculations were performed on the three models using regular grids of 256x256 square elements. 
Table 6 shows the effective Young’s modulus and Poisson’s ratio obtained from the three RVEs. The FEM estimates 
of Young’s modulus are 9 times higher than the Young’s modulus of PBX 9501. The higher stiffness obtained from 
the FEM calculations is partly because the large number of contacts between particles lead essentially to a continuous 
particle phase containing pockets of binder. 

Micrographs of PBX 9501 show that there is considerable damage in the particles, binder and the particle-binder 
interface [16], Better estimates of the effective properties can probably be obtained by incorporating damage in the 
models of PBX 9501 discussed above. However, this issue is not explored in this work because the amount and type 
of damage is difficult to quantify in polymer bonded explosives. 

5 Effective moduli of PBX 9501 

Of the various microstructures of PBX 9501 explored in the previous section, the best estimate of the Young’s modulus 
of PBX 9501 is provided by the 100 particle model of the dry blend of PBX 9501 shown in Figure 8, overlaid 
by a 350x350 grid. The effective properties of PBX 9501 shown in this section have been calculated at different 
temperatures and strain rates using the 100 particle model of Figure 8. Instead of the experimentally determined 
values of moduli of HMX, data from molecular dynamics (MD) simulations [18] have been used because estimates 
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using MD data of the Young’s modulus of PBX 9501 at room temperature and low strain rate are slightly closer to the 
experimental value than estimates calculated using the experimentally determined modulus of HMX. 

The Young’s modulus of HMX from MD simulations is 17,700 MPa and the Poisson’s ratio is 0.21 compared to a 
Young’s modulus of 15,300 MPa and a Poisson’s ratio of 0.32 from experimental data. The Young’s modulus of HMX 
is assumed constant at the temperatures and strain rates considered. 

The PBX 9501 binder data collected from various sources [19, 20, 6, 21] are listed in Table 7. The data are sorted 
in order of decreasing temperature followed by increasing strain rate. As shown in Table 7, the Young’s modulus of the 
binder generally increases with increase in strain rate and with decrease in temperature though some of the available 
data do not follow this trend. Since values for the Poisson’s ratio of the binder are not available, a value of 0.49 has 
been used for all cases. However, the value is likely to be lower near the glass transition temperature of estane (around 
-33°C). As before, a dirty binder has been used for the grid cells not fully occupied by particles. The differential 
effective medium approximation [2] has been used to calculate the effective properties of the dirty binder using the 
moduli of HMX and the binder. Table 7 also shows the dirty binder properties used in the simulations. It is interesting 
to note that the Poisson’s ratio of the dirty binder does not change significantly even though the Young’s modulus can 
vary considerably depending on the modulus of the PBX 9501 binder. 

Table 8 shows the effective properties of the model RVE computed using finite elements along with the available 
experimental data on PBX 9501 . The experimental data are within 15% of the finite element estimates for temperatures 
close to 25°C and at low strain rates. Good agreement is also observed for the data near 16°C and at high strain rates. 
However, at 0°C and a strain rate of 2200/s the experimental value is around 2 times the predicted values of Young’s 
modulus while at -20°C and high strain rate, the experimental value is half the estimate. Experimental data on the 
Poisson’s ratio of PBX 9501 are available only at 25°C and 0.005/s strain rate. The predicted Poisson’s ratio under 
these conditions is considerably lower that the experimental value. Since the Poisson’s ratio is highly sensitive to 
numerical error [22], more accurate numerical calculations may be required to resolve this issue. 

These results show that, based on existing data on PBX 9501, it is not possible to conclude that the effective moduli 
of PBX 9501 can be predicted accurately if the moduli of its constituents are known for various temperatures and strain 
rates. However, the results do show that for certain RVEs, close approximations to the effective Young’s modulus of 
PBX 9501 can be obtained from two-dimensional finite element simulations. 
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6 Summary and conclusions 



A two-dimensional finite element based method for approximating the effective elastic moduli of particulate compos- 
ites has been presented. This approach has been applied to representative volume elements (RVEs) containing circular 
particles for volume fractions from 0.1 to 0.92. The estimated Young’s moduli and Poisson’s ratios have been found to 
closely approximate those from differential effective medium approximations. Comparisons with three-dimensional 
finite element simulations have also shown that the two-dimensional approach provides good approximations for RVEs 
containing polydisperse spherical particles. The estimates of Young’s modulus from three-dimensional models have 
been found to be slightly higher than from two-dimensional models. On the other hand, the Poisson’s ratios from 
three-dimensional models have been found to be lower than estimates from two-dimensional models. 

Since polymer bonded explosives such as PBX 9501 contain two widely separated particle sizes, a number of 
manually generated microstructures of PBX 9501 have been simulated. The models were designed to contain a few 
large particles and a larger number of smaller interstitial particles while avoiding particle -particle contact. Such models 
have been found to underestimate the Young’s modulus of PBX 9501 considerably when discretized with triangular 
elements. However, when the same models are discretized with square elements, considerably different estimates of 
Young’s modulus are obtained depending on the degree of discretization. Hence, it is extremely difficult to determine 
an appropriate RVE for high modulus contrast and high volume fraction particulate composites that has the optimal 
size, number of particles, particle shapes and size distributions. The RVE also has to be such that it can be easily 
discretized using elements that are accurate and relatively computationally inexpensive. 

Randomly generated microstructures based on the particle size distribution of PBX 9501 show that the size of 
the RVE plays a lesser role in the determination of the effective Young’s modulus than the amount of discretization 
used, even for highly refined meshes. If aligned square particles are used instead of circular particles, the predicted 
moduli are considerably higher than the experimental moduli because contacts between particles increase dramatically. 
It is also observed that the numerical simulations consistently underestimate the effective Poisson’s ratio. Since the 
Poisson’s ratio is highly sensitive to numerical error, accurate numerical simulations have to be performed to resolve 
this issue. 

Simulations of PBX 9501 for various temperatures and strain rates suggest that if a RVE is chosen that predicts 
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a reasonable value of Young’s modulus at a given temperature and strain rate, it can be used to obtain acceptable 
estimates at other temperatures and strain rates. Further experimental data are required to confirm this finding. 
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A Two- and Three-Dimensional Moduli 

The determination of the effective properties of random particulate composites using numerical techniques requires 
the use of a three-dimensional representative volume element (RVE). If a two-dimensional RVE is selected, any plane 
strain representation immediately implies that the material is composed of unidirectional cylindrical particles. On the 
other hand, if a plane stress representation is used, the physical interpretation of the composite is a thin sheet cut from 
the particulate composite. 

One of the intents of this work has been to show that acceptable estimates of the effective elastic moduli of par- 
ticulate composites can be obtained from two-dimensional plane strain calculations of the stress and strain fields. 
However, two-dimensional calculations in a single plane of the composite lead to two-dimensional or planar effective 
elastic moduli. For example, if a uniform uniaxial displacement Ui is applied to a unit square of an isotropic, homo- 
geneous, linear elastic material with Young’s modulus E and Poisson’s ratio v in plane strain, the resulting reaction 
force is F\ and the displacement perpendicular to the direction of the applied displacement is — Since the square 
is of unit size, the nonzero strains are qi = u i and £22 = — 1 * 2 - The nonzero stresses are an = F t and a 33 = vF\ , 
the latter obtained from the constitutive relation 

= (1/E)[(1 + v)a ij - vo kk 5, tj \ i,j,k = 1,2,3. (6) 

Note that E and v are three-dimensional or true linear elastic moduli. It can be shown that the apparent Poisson’s ratio 
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in the plane (€22/^1 1) is given by 



V 2 D = £22/611 = u 2 /ui = v/(l - v) (7) 

and the apparent Young’s modulus in the plane (cth/ch) is given by 

E- 2 D = 0'11/ell = Fi/ui = E/( 1 — v 2 ) (8) 

These plane moduli, E 2 D and v 2 u • are referred to as two-dimensional moduli following Jun and Jasiuk [3] and 
Torquato [23], The reason for such a name is that these apparent moduli are the same as those obtained if the consti- 
tutive equations are expressed, in two-dimensional form, as 

Cij = (1/E 2D )[(1 + v 2 D)&ij - v 2 D^kAj\ i,j, k = 1, 2. (9) 

The results presented in this work suggest that these relations between two- and three-dimensional elastic mod- 
uli (which are exact only for isotropic, homogeneous, linear elastic materials) can be used to obtain reasonable ap- 
proximations of three-dimensional elastic moduli of random particulate composites from two-dimensional numerical 
simulations. 
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Figure captions 



Figure 1. RVEs containing 10% to 92% by volume of circular particles. f p is the volume fraction of particles in a 
RVE. 

Figure 2. Comparison of finite element (FEM) and differential effective medium (DEM) predictions. E is the Young’s 
modulus of the binder. 

Figure 3. Two- and three-dimensional RVEs containing 70%, 75% and 80% particles by volume. f p is the volume 
fraction of particles. 

Figure 4. Two- and three-dimensional estimates of effective Young’s modulus and Poisson’s ratio. 

Figure 5. Microstructure of PBX 9501 (adapted from [7]). The largest particles are around 300 microns in size. 
Figure 6. Manually generated microstructures containing ~ 90% by volume of circular particles. 

Figure 7. Microstructure containing 92% particles modeled with six-noded triangles and four-noded squares. 

Figure 8. Circular particle models based on the dry blend of PBX 9501 . 

Figure 9. Circular particle models based on the particle size distribution of pressed PBX 9501. 

Figure 10. Square particle model based on the size distribution of pressed PBX 9501. 
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Tables 
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Table 1: Effective moduli of the six model PBX 9501 microstructures from FEM calculations six-noded triangular 
elements. E is the Young’s modulus and v is the Poisson’s ratio. 





Expt. 


1 


2 


Model RVE 
3 4 


5 


6 


Mean 


Std. Dev. 


E (MPa) 


1013 


116 


126 


130 


42 


183 


192 


132 


54 


V 


0.35 


0.34 


0.32 


0.32 


0.44 


0.28 


0.25 


0.33 


0.07 
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Table 2: Effective elastic moduli of the four models of the dry blend of PBX 9501. 



No. of Size Young’s Modulus (MPa) Poisson’s Ratio 

Particles (mm) FEM Expt. FEM Expt. 

256x256 350x350 256x256 350x350 



100 


0.65 


1959 


968 




0.22 


0.20 


0.35 






2316 


1488 




0.23 


0.23 


0.35 






2899 


2004 




0.25 


0.24 


0.35 






4350 


2845 




0.25 


0.25 


0.35 
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Table 3: Material properties of the dirty binder for the four pressed PBX 9501 microstructures. f p is the volume 
fraction of particles. 



No. of 
Particles 


fp 


fp 

in binder 


Properties of dirty binder 
E (MPa) v 


100 


0.89 


0.28 


1.583 


0.484 


200 


0.87 


0.37 


2.1395 


0.481 


500 


0.86 


0.43 


2.713 


0.478 


1000 


0.855 


0.45 


2.952 


0.477 
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Table 4: Effective elastic moduli of the four models of pressed PBX 9501. 



No. of Size Young’s Modulus (MPa) Poisson’s Ratio 

Particles (mm) FEM Expt. FEM Expt. 

256x256 350x350 256x256 350x350 



100 


0.36 


2925 


1798 




0.23 


0.14 


0.35 




0.42 


3342 


2408 




0.21 


0.20 


0.35 


500 


0.54 


5256 


3994 




0.25 


0.23 


0.35 


1000 


0.68 


6171 


4756 




0.25 


0.25 


0.35 
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Table 5: Material properties of the dirty binder for the three PBX microstructures with square particles. 



No. of f p f p Properties of dirty binder 

Particles in binder E (MPa) v 



700 




0.39 


2.358 


0.4799 


2800 




0.40 


2.448 


0.4795 


11600 




0.42 


2.588 


0.4788 
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Table 6: Effective stiffnesses of the model microstructures with square particles. 



No. of Size Young’s Modulus (MPa) Poisson’s Ratio 

Particles (mm) FEM (256 x 256) Expt. FEM (256 x 256) Expt. 



700 


3.6 


9119 




0.26 


0.35 


2800 


5.3 


9071 


1013 


0.27 


0.35 


11600 


9.0 


9593 


1013 


0.27 


0.35 
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Table 7: Material properties used to calculate the effective moduli of PBX 9501 for different temperatures and strain 
rates. Numbers in square brackets refer to the source of the binder data. 



Temp. 

(°C) 


Strain 

Rate 

(/s) 


HMX 

Young’s Poisson’s 
Modulus Ratio 

(MPa) 


PBX 9501 Binder 
Young’s Poisson’s 
Modulus Ratio 

(MPa) 


’Dirty’ Binder 
Young’s Poisson’s 
Modulus Ratio 

(MPa) 


25 


0.005 


17700 


0.21 


0.59 [6] 


0.49 


2.15 


0.479 


25 


0.008 


17700 


0.21 


0.73 [6] 


0.49 


2.66 


0.479 


25 


0.034 


17700 


0.21 


0.81 [6] 


0.49 


2.95 


0.479 


25 


0.049 


17700 


0.21 


0.82 [6] 


0.49 


2.99 


0.479 


25 


2400 


17700 


0.21 


300 [20] 


0.49 


1001 


0.474 


22 


0.001 


17700 


0.21 


0.47 [21] 


0.49 


1.71 


0.479 


22 


1 


17700 


0.21 


1.4 [21] 


0.49 


5.10 


0.479 


22 


2200 


17700 


0.21 


3.3 [21] 


0.49 


12.0 


0.479 


16 


1700 


17700 


0.21 


22.5 [19] 


0.49 


81.5 


0.479 


0 


0.001 


17700 


0.21 


0.85 [21] 


0.49 


3.10 


0.479 


0 


1700 


17700 


0.21 


246 [19] 


0.49 


833 


0.475 


0 


2200 


17700 


0.21 


4 [21] 


0.49 


14.6 


0.479 


-15 


0.001 


17700 


0.21 


1.4 [21] 


0.49 


5.10 


0.479 


-15 


1 


17700 


0.21 


5.7 [21] 


0.49 


20.7 


0.479 


-15 


1000 


17700 


0.21 


1600 [21] 


0.49 


4082 


0.458 


-20 


0.001 


17700 


0.21 


1.6 [21] 


0.49 


5.83 


0.479 


-20 


1200 


17700 


0.21 


1600 [21] 


0.49 


4082 


0.458 


-20 


1700 


17700 


0.21 


1333 [19] 


0.49 


3556 


0.461 


-40 


0.001 


17700 


0.21 


5.7 [21] 


0.49 


20.7 


0.479 


-40 


0.001 


17700 


0.21 


5.3 [19] 


0.49 


19.3 


0.479 


-40 


1300 


17700 


0.21 


1000 [21] 


0.49 


2838 


0.464 
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Table 8: Effective Young’s modulus from finite element simulations of a 100 particle model of the dry blend of PBX 
9501. Numbers in square brackets refer to the source of the PBX 9501 data. 



Temp. 

(°C) 


Strain 

Rate 

(/s) 


FEM Estimates 
Young’s Poisson’s 
Modulus Ratio 

(MPa) 


Experiments 
Young’s Poisson’s 

Modulus Ratio 

(MPa) 


25 


0.005 


1001 


0.19 


1040 [24] a 


- 


25 


0.01 


1080 


0.20 


770 [20] 


- 










1020 [19] b 


- 


25 


0.034 


1123 


0.20 


- 


- 


25 


0.05 


1129 


0.20 


1013 [6] 


0.35 


25 


2400 


11273 


0.31 


- 


- 


22 


0.001 


926 


0.17 


- 


- 


22 


1 


1399 


0.23 


- 


- 


22 


2200 


2050 


0.26 


- 


- 


16 


1700 


4981 


0.30 


4650 [19] c 


- 


0 


0.001 


1144 


0.21 


- 


- 


0 


1700 


10808 


0.31 


- 


- 


0 


2200 


2243 


0.27 


5480 [19] 


- 


-15 


0.001 


1399 


0.23 


- 


- 


-15 


1 


2652 


0.28 


- 


- 


-15 


1000 


14602 


0.28 


- 


- 


-20 


0.001 


1481 


0.23 


- 


- 


-20 


1200 


14602 


0.28 


- 


- 


-20 


1700 


14291 


0.29 


6670 [19] d 


- 


-40 


0.001 


2652 


0.28 


- 


- 


-40 


0.001 


2562 


0.27 


- 


- 


-40 


1300 


13778 


0.29 


12900 [19] d 


- 



a - experimental data are for a strain rate of 0.001/s. 

b - experimental data are for a strain rate of 0.01 1/s and a temperature of 27°C. 
c - experimental data are for a strain rate of 2250/s and a temperature of 17°C. 
d - experimental data are for a strain rate of 2250/s. 
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f =0.7 f =0.8 f =0.92 

p p p 



Figure 1: RVEs containing 10% to 92% by volume of circular particles. f p is the volume fraction of particles in a 
RVE. 
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Figure 1. 

Biswajit Banerjee and Daniel O. Adams 



34 





0.5j 

1 






0.45 








N 


V \\ 


0.4 




O FEM(E b =l) 


\v 






DEM(E b = 1) 








x FEM (E = 10) 


\ V '' \ 


0.35 




DEM (E fc = 10) 


\ 'A 






n FEM (E = 100) 


\ A 


0.3 




DEM (E =100) 


\ A ■ 






FEM (E = 1000) 








DEM (E fe = 1000) 


\h 


0.25 




FEM (E b = 10000) 


Ya 






DEM = 10000) 


A 


0.2 


L 




L i \t 



0.2 u L- ■- 

0 0.2 0.4 0.6 0.8 

Particle Vol. Frac. 



(b)Poisson’s Ratio. 



Figure 2: Comparison of finite element (FEM) and differential effective medium (DEM) predictions. E(, is the Young’s 
modulus of the binder. 
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Figure 2. 

Biswajit Banerjee and Daniel O. Adams 
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fp = 0.7 



f P = 0.75 



f P = 0.8 





(a) Two-dimensional RVEs. 



f P = 0.7 f p = 0.75 f p = 0.8 




(b) Three-dimensional RVES. 



Figure 3: Two- and three-dimensional RVEs containing 70%, 75% and 80% particles by volume. f p is the volume 
fraction of particles. 
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Figure 3. 

Biswajit Banerjee and Daniel O. Adams 
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Young's Modulus (MPa) 




Particle Vol. Frac. Particle Vol. Frac. 

(a) Young’s Modulus. (b) Poisson’s Ratio. 

Figure 4: Two- and three-dimensional estimates of effective Young’s modulus and Poisson’s ratio 






Figure 4. 

Biswajit Banerjee and Daniel O. Adams 



40 





Figure 5: Microstructure of PBX 9501 (adapted from [7]). The largest particles are around 300 microns in size. 
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Figure 5. 

Biswajit Banerjee and Daniel O. Adams 
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Model 4 Model 5 Model 6 

Figure 6: Manually generated microstructures containing ~ 90% by volume of circular particles. 
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Figure 6. 

Biswajit Banerjee and Daniel O. Adams 
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(c) Six-noded Triangles. (d) Four-noded Squares. 



Figure 7: Microstructure containing 92% particles modeled with six-noded triangles and four-noded squares. 
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Figure 7. 

Biswajit Banerjee and Daniel O. Adams 
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0.65x0.65 mm 2 
100 Particles 





0.94x0.94 mm 2 1.13x1.13mm 2 

200 Particles 300 Particles 




1.33x 1.33 mm 2 
400 Particles 



Figure 8: Circular particle models based on the dry blend of PBX 9501. 
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Figure 8. 

Biswajit Banerjee and Daniel O. Adams 
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0.36x0.36 mm 2 
100 Particles 




0.42x0.42 mm 2 0.54x0.54 mm 2 

200 Particles 500 Particles 




0.68x0.68 mm 2 
1000 Particles 



Figure 9: Circular particle models based on the particle size distribution of pressed PBX 9501. 
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Figure 9. 

Biswajit Banerjee and Daniel O. Adams 
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9.0x9.0 mm 2 
11600 Particles 



Figure 10: Square particle model based on the size distribution of pressed PBX 9501. 



3. 6x3. 6 mm 2 5. 3x5. 3 mm 2 

700 Particles 2800 Particles 
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Figure 10. 

Biswajit Banerjee and Daniel O. Adams 
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